
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The impact of divorce laws on the equilibrium in the marriage market.
% Ana Reynoso
% April 2024
%
% This file inputs estimated parameters and moments under mutual consent
% regime and produces:
%
% - Figure 3,
% - Table 1,
% - Table A4,
% - Table 2,
% - Figure 4, and
% - Figure A1.
%  
% Data: PSID 1968-1992
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%------------------------ Replication Path -------------------------------%

clc
clear all

%--- Indicate location of Replication folder:
%replication_location = 'C:\update_with_your_path';

%------------------------ Preliminaries ----------------------------------%

%--- Working directory
cd ([replication_location,'\Figures_and_Tables']);

%--- Path for data
addpath(genpath([replication_location,'\Data\Inputs']));

%--- Path for functions
addpath(genpath([replication_location,'\Figures_and_Tables\matrix2latex']));

%--- Path for outputs
output_dir = [replication_location,'\Outputs'];

%--- Parameters for figures
x0=0;
y0=0;
width=1800;
height=1200;

%--- Load estimation outputs
load([output_dir, '\MCD\outputs_mcd']);
params=X;
load([output_dir, '\MCD\Variance_mcd']);
load([output_dir, '\MCD\Sensitivity_mcd']);
load([output_dir, '\MCD\moments_mc']);
load([output_dir, '\MCD\moments_mc_untargeted']);
                                   
%------------------------ Figure 3 ---------------------------------------%

female_educ_bymkt = [data_f(1,:); data_f(4,:); data_f(7,:)];
male_educ_bymkt   =  data_m(1:3,:);

pop_m1 = [female_educ_bymkt(:,1) male_educ_bymkt(:,1)];
pop_m2 = [female_educ_bymkt(:,2) male_educ_bymkt(:,2)];
pop_m3 = [female_educ_bymkt(:,3) male_educ_bymkt(:,3)];
pop_m4 = [female_educ_bymkt(:,4) male_educ_bymkt(:,4)];

female_shr_bymkt = round([female_educ_bymkt(2,:)./ female_educ_bymkt(1,:); female_educ_bymkt(3,:)./ female_educ_bymkt(1,:)],2);
male_shr_bymkt = round([male_educ_bymkt(2,:)./ male_educ_bymkt(1,:); male_educ_bymkt(3,:)./ male_educ_bymkt(1,:)],2);

tiledlayout(2,2);

%--- Mkt1
nexttile;
name={'hs';'sc';'c+'};

p1str=num2str(female_shr_bymkt(1,1));
p2str=num2str(male_shr_bymkt(1,1)) ;
p3str=num2str(female_shr_bymkt(2,1)) ;
p4str=num2str(male_shr_bymkt(2,1)) ;

labels1 = {p1str, p3str};
labels2 = {p2str, p4str};

a = bar(pop_m1, 'grouped');
xlabel('Northeast')
legend('Women','Men','Location','northwest','FontSize',20, 'NumColumns', 2) 
set(gca,'xticklabel',name,'FontSize',30, 'FontWeight','bold');
ylim([0 0.65])
a(1).FaceColor = [0 0.3470 0.5410]; 
a(2).FaceColor = [0.9 0.9 0.9]; 
text([1.95 2.95], pop_m1(2:3,1), labels1, 'HorizontalAlignment','right', 'VerticalAlignment','bottom','FontSize',30, 'FontWeight','bold')
text([2.05 3.05], pop_m1(2:3,2), labels2, 'HorizontalAlignment','left', 'VerticalAlignment','bottom','FontSize',30, 'FontWeight','bold')

%--- Mkt2
nexttile;

p1str=num2str(female_shr_bymkt(1,2)) ;
p2str=num2str(male_shr_bymkt(1,2)) ;
p3str=num2str(female_shr_bymkt(2,2));
p4str=num2str(male_shr_bymkt(2,2)) ;

labels1 = {p1str, p3str};
labels2 = {p2str, p4str};

a = bar(pop_m2, 'grouped');
xlabel('Midwest and West')
set(gca,'xticklabel',name,'FontSize',30, 'FontWeight','bold');
ylim([0 0.65])
a(1).FaceColor = [0 0.3470 0.5410];
a(2).FaceColor = [0.9 0.9 0.9];
text([1.95 2.95], pop_m2(2:3,1), labels1, 'HorizontalAlignment','right', 'VerticalAlignment','bottom','FontSize',30, 'FontWeight','bold')
text([2.05 3.05], pop_m2(2:3,2), labels2, 'HorizontalAlignment','left', 'VerticalAlignment','bottom','FontSize',30, 'FontWeight','bold')

%--- Mkt3
nexttile;

p1str=num2str(female_shr_bymkt(1,3));
p2str=num2str(male_shr_bymkt(1,3));
p3str=num2str(female_shr_bymkt(2,3));
p4str=num2str(male_shr_bymkt(2,3));

labels1 = {p1str, p3str};
labels2 = {p2str, p4str};

a = bar(pop_m3, 'grouped');
xlabel('South Atlantic')
set(gca,'xticklabel',name,'FontSize',30, 'FontWeight','bold');
ylim([0 0.65])
a(1).FaceColor = [0 0.3470 0.5410];
a(2).FaceColor = [0.9 0.9 0.9];
text([1.95 2.95], pop_m3(2:3,1), labels1, 'HorizontalAlignment','right', 'VerticalAlignment','bottom','FontSize',30, 'FontWeight','bold')
text([2.05 3.05], pop_m3(2:3,2), labels2, 'HorizontalAlignment','left', 'VerticalAlignment','bottom','FontSize',30, 'FontWeight','bold')

%--- Mkt4
nexttile;

p1str=num2str(female_shr_bymkt(1,4)) ;
p2str=num2str(male_shr_bymkt(1,4)) ;
p3str=num2str(female_shr_bymkt(2,4));
p4str=num2str(male_shr_bymkt(2,4));

labels1 = {p1str, p3str};
labels2 = {p2str, p4str};

a = bar(pop_m4, 'grouped');
xlabel('South Central')
set(gca,'xticklabel',name,'FontSize',30, 'FontWeight','bold');
ylim([0 0.65])
a(1).FaceColor = [0 0.3470 0.5410];
a(2).FaceColor = [0.9 0.9 0.9];
text([1.95 2.95], pop_m4(2:3,1), labels1, 'HorizontalAlignment','right', 'VerticalAlignment','bottom','FontSize',30, 'FontWeight','bold')
text([2.05 3.05], pop_m4(2:3,2), labels2, 'HorizontalAlignment','left', 'VerticalAlignment','bottom','FontSize',30, 'FontWeight','bold')

%--- Store
set(gcf,'position',[x0,y0,width,height])
saveas(gcf, 'Output/Figure_3.png')  
close

%------------------------ Table 1 ----------------------------------------%

%--- Sample size
data_imp = importdata('sample_sizes_4mkts.csv');
data_sampsize = data_imp.data;

samp_size_m1 = data_sampsize(1:15,:);
samp_size_m2 = data_sampsize(16:30,:);
samp_size_m3 = data_sampsize(31:45,:);
samp_size_m4 = data_sampsize(46:60,:);

samp_size_all = [samp_size_m1(:,2) samp_size_m1(:,3)+samp_size_m2(:,3)+samp_size_m3(:,3)+samp_size_m4(:,3) samp_size_m1(:,4)+samp_size_m2(:,4)+samp_size_m3(:,4)+samp_size_m4(:,4) ];
samp_size_total = sum(samp_size_all);

%--- Arrange estimates, estandard errors, and sensitivity moments
param_ests=zeros(9,3); 
param_ests(:,1)=params(1:9).*[params(37) params(37) params(37) params(38) params(38) params(38) params(39) params(39) params(39)]; 
param_ests(:,2)=params(10:18); 
single_ests=zeros(3,2); 
single_ests(:,1)=params(19:21);
single_ests(:,2)=params(22:24);
param_ests(:,3)=params(25:33); 

se_couples=zeros(15,3);
se_couples(1:9,1)=se(1:9);
se_couples(1:9,2)=se(10:18);
se_couples(10:15,2)=se(19:24);
se_couples(1:9,3)=se(25:33);

est_se=zeros(15,11);
est_se(:,1)= samp_size_all(:,2);
est_se(:,2)= samp_size_all(:,3);
est_se(1:9,3)= param_ests(:,1);
est_se(1:9,6)= param_ests(:,2);
est_se(10:12,6)= single_ests(:,1);
est_se(13:15,6)= single_ests(:,2);
est_se(1:9,9)= param_ests(:,3);
est_se(:,4)= se_couples(:,1);
est_se(:,10)= se_couples(:,3);
est_se(:,7)= se_couples(:,2);
est_se = round(est_se,2);

mom_sahw = cell(15,3);
mom_sahw(1:9,:)=sensitivity_moments(1:9,:);
mom_sahw(10:end,:)={'nan'};
mom_mq = sensitivity_moments(10:24,:);
mom_var = cell(15,3);
mom_var(1:9,:)=sensitivity_moments(25:33,:);
mom_var(10:end,:)={'nan'};

%--- Build the table
est_se_mom = table(est_se(:,1),est_se(:,2),est_se(:,3),est_se(:,4),mom_sahw,...                                   
   est_se(:,6),est_se(:,7), mom_mq, est_se(:,9),est_se(:,10),mom_var, 'VariableNames', {'Hhs.', 'Obs.', 'S-a-h Est.','S-a-h s.e', 'S-a-h Moments',  'Mean match quality Est.', ...
   'Mean match quality s.e.','Mean match quality Moments', 'Variance match quality Est.', 'Variance match quality s.e', 'Variance match quality Moments'}, ...
    'RowNames', {'(hs,hs)', '(hs,sc)', '(hs,c$+$)','(sc,hs)', '(sc,sc)', '(sc,c$+$)','(c$+$,hs)','(c$+$,sc)', '(c$+$,c$+$)','hs (females)','sc (females)','c$+$ (females)','hs (males)','sc (males)','c$+$ (males)'});

%--- Store
writetable(est_se_mom, 'Output/Table_1.xls')

%------------------------ Table A4 ---------------------------------------%

%--- Arrange observations, estimates, estandard errors, and sensitivity moments
param_ests_w=zeros(6,4);
param_ests_w(:,1)=params(34:39); 
param_ests_w(:,2)=params(40:45); 
param_ests_w(:,3)=params(46:51); 
param_ests_w(1:3,4)=params(52:54); 

se_wages=zeros(6,4);
se_wages(:,1)=se(34:39); 
se_wages(:,2)=se(40:45); 
se_wages(:,3)=se(46:51); 
se_wages(1:3,4)=se(52:54); 

est_se_w= [param_ests_w(:,1) se_wages(:,1) param_ests_w(:,2) se_wages(:,2) param_ests_w(:,3) se_wages(:,3) param_ests_w(:,4) se_wages(:,4) ];
est_se_w(:,[1 2]) = round(est_se_w(:,[1 2]),2);
est_se_w(:,3:end) = round(est_se_w(:,3:end),4);

mom_price = cell(6,3);
mom_price(1:6,:)=sensitivity_moments(34:39,:);
mom_exp = cell(6,3);
mom_exp(1:6,:)=sensitivity_moments(40:45,:);
mom_exp2 = cell(6,3);
mom_exp2(1:6,:)=sensitivity_moments(46:51,:);
mom_sp = cell(6,3);
mom_sp(1:3,:)=sensitivity_moments(52:54,:);
mom_sp(4:end,:)={'nan'};

obs_w= [ est_se(1,2) + est_se(4,2) + est_se(7,2) + est_se(13,2) ; ...
         est_se(2,2) + est_se(5,2) + est_se(8,2) + est_se(14,2) ; ...
         est_se(3,2) + est_se(6,2) + est_se(9,2) + est_se(15,2) ; ...
         est_se(1,2) + est_se(2,2) + est_se(3,2) + est_se(10,2) ; ...
         est_se(4,2) + est_se(5,2) + est_se(6,2) + est_se(11,2) ; ...
         est_se(7,2) + est_se(8,2) + est_se(9,2) + est_se(12,2) ]; 

%--- Build the table
est_se_mom_w = table(obs_w, est_se_w(:,1), est_se_w(:,2),mom_price,...                       %mean price
   est_se_w(:,3),est_se_w(:,4),mom_exp, est_se_w(:,5),est_se_w(:,6),mom_exp2,...     %exp and exp2
   est_se_w(:,7),est_se_w(:,8),mom_sp, ... %exp and exp2
   'VariableNames', {'Obs.', 'Mean price of education Est.','Mean price of education s.e', 'Mean price of education Moments',  'experience Est.', 'experience s.e.', ...
   'experience Moments', 'experience2 Est.', 'experience2 s.e', 'experience2 Moments', 'spousal experience Est.', 'spousal experience s.e', 'spousal experience Moments'}, ...
    'RowNames', {'hs (men)','sc (men)','c$+$ (men)','hs (women)','sc (women)','c$+$ (women)'})   ;

%--- Store
writetable(est_se_mom_w, 'Output/Table_A4.xls')

%------------------------ Table 2 ----------------------------------------%

lambda_mc = reshape(params(55:90), 9, 4);

pweight_mcd_avg = round(lambda_mc(:,1)*0.23 + lambda_mc(:,2)*0.23+ lambda_mc(:,3)*0.36 + lambda_mc(:,4)*0.18,2);
pweight_mc_avg_table=[pweight_mcd_avg(1:3)';pweight_mcd_avg(4:6)';pweight_mcd_avg(7:9)'];
rowLabels = {'hs','sc','c$+$'};
columnLabels = {'hs','sc','c$+$'};
matrix2latex(pweight_mc_avg_table, 'Output/Table_2.tex',  'rowLabels', rowLabels, 'columnLabels', columnLabels, 'alignment', 'c', 'format', '%-6.2f', 'size', 'tiny');

%------------------------ Figure 4 ---------------------------------------%

%--- Model and data moments
model_moments_bymkt_mc = reshape(model_moments_mc(1:168), 168/4,4); 
data_moments_bymkt = reshape(data_moments(1:168), size(model_moments_bymkt_mc,1),size(model_moments_bymkt_mc,2)); 

%--- Confidence intervals
mom_ci = importdata('moments_boot_ci.csv');
ci_mf = mom_ci.data(:,1:18);
ci_mf_l = ci_mf(:,[1 3 5 7 9 11 13 15 17]);
ci_mf_u = ci_mf(:,[2 4 6 8 10 12 14 16 18]);

%--- Build figure
figure(1)
set(gcf,'Color','White')
b=bar(model_moments_bymkt_mc(1:9,:), 'grouped');
hold on

[ngroups,nbars] = size(model_moments_bymkt_mc(1:9,:));
xpos = nan(nbars, ngroups);
for i = 1:nbars
    xpos(i,:) = b(i).XEndPoints;
end

%mkt 1
plot(xpos(1,:),data_moments_bymkt(1:9,1)','o','MarkerSize',8,'MarkerEdgeColor','red','MarkerFaceColor','red');
errorbar( xpos(1,:)',  data_moments_bymkt(1:9,1)',                ...
            data_moments_bymkt(1:9,1)'-ci_mf_l(1,:),                                    ...
            ci_mf_u(1,:)-data_moments_bymkt(1:9,1)',                                    ...
            'vertical','.','MarkerSize',0.01,'LineWidth',2,'Color','red' );
hold on       
%mkt 2
plot(xpos(2,:),data_moments_bymkt(1:9,2)','o','MarkerSize',8,'MarkerEdgeColor','red','MarkerFaceColor','red');
errorbar( xpos(2,:)',  data_moments_bymkt(1:9,2)',                ...
            data_moments_bymkt(1:9,2)'-ci_mf_l(2,:),                                    ...
            ci_mf_u(2,:)-data_moments_bymkt(1:9,2)',                                    ...
            'vertical','.','MarkerSize',0.01,'LineWidth',2,'Color','red' );
hold on       
%mkt 3
plot(xpos(3,:),data_moments_bymkt(1:9,3)','o','MarkerSize',8,'MarkerEdgeColor','red','MarkerFaceColor','red');
errorbar( xpos(3,:)',  data_moments_bymkt(1:9,3)',                ...
            data_moments_bymkt(1:9,3)'-ci_mf_l(3,:),                                    ...
            ci_mf_u(3,:)-data_moments_bymkt(1:9,3)',                                    ...
            'vertical','.','MarkerSize',0.01,'LineWidth',2,'Color','red' );
hold on       
%mkt 4
plot(xpos(4,:),data_moments_bymkt(1:9,4)','o','MarkerSize',8,'MarkerEdgeColor','red','MarkerFaceColor','red');
errorbar( xpos(4,:)',  data_moments_bymkt(1:9,4)',                ...
            data_moments_bymkt(1:9,4)'-ci_mf_l(4,:),                                    ...
            ci_mf_u(4,:)-data_moments_bymkt(1:9,4)',                                    ...
            'vertical','.','MarkerSize',0.01,'LineWidth',2,'Color','red' );
       
xticks(1:9);
xticklabels({   '(hs,hs)','(hs,sc)',        ...
                '(hs,c+)','(sc,hs)','(sc,sc)', '(sc,c+)', '(c+,hs)', ...
                '(c+,sc)','(c+,c+)'})
ax                          = gca;
ax.FontSize                 = 25;
ax.FontWeight               ='bold';
b(1).FaceColor = [0 0 0]; 
b(2).FaceColor = [0 0.3470 0.5410]; 
b(3).FaceColor = [0.9 0.9 0.9]; 
b(4).FaceColor = [1 1 1];

legend('Northeast','Midwest and West','South Atlantic','South Central','Data','Confidence Interval','Location','northeast','FontSize',20, 'FontWeight','bold','NumColumns',3) ;

%--- Store
set(gcf,'position',[x0,y0,width,height]);
saveas(gcf, 'Output/Figure_4.png')  ;
close

%------------------------ Figure A1 --------------------------------------%

%--- Model and data moments
data_imp_utg = importdata('data_untargeted.csv');
data_utg = data_imp_utg.data;
                                
duration_all_data =data_utg(:,2);
data_sahw_t = data_utg(:,3:end);

data_moments_bymkt = reshape(data_moments(1:168), 168/4,4); 
women_moments_bymkt = data_moments_bymkt(1:9,:)./data_f;

data_sahw_t_byeduct = [sum(data_sahw_t(1:3,:,:).*women_moments_bymkt(1:3,1)); ...
                       sum(data_sahw_t(4:6,:,:).*women_moments_bymkt(4:6,1));...
                      sum(data_sahw_t(7:9,:,:).*women_moments_bymkt(7:9,1))];

%--- Confidence intervals
data_imp_utg_CI = importdata('untg_moments_boot_CI.csv');
data_utg_CI = data_imp_utg_CI.data;

duration_lb = data_utg_CI(2,1:9);
duration_ub = data_utg_CI(1,1:9);

data_sahw_t_lb = zeros(9,4);
data_sahw_t_lb(:,1) = data_utg_CI(2,10:18)';
data_sahw_t_lb(:,2) = data_utg_CI(2,19:27)';
data_sahw_t_lb(:,3) = data_utg_CI(2,28:36)';
data_sahw_t_lb(:,4) = data_utg_CI(2,37:45)';

data_sahw_t_ub = zeros(9,4);
data_sahw_t_ub(:,1) = data_utg_CI(1,10:18)';
data_sahw_t_ub(:,2) = data_utg_CI(1,19:27)';
data_sahw_t_ub(:,3) = data_utg_CI(1,28:36)';
data_sahw_t_ub(:,4) = data_utg_CI(1,37:45)';

data_sahw_t_byeduct_lb = [sum(data_sahw_t_lb(1:3,:,:).*women_moments_bymkt(1:3,1)); ...
                       sum(data_sahw_t_lb(4:6,:,:).*women_moments_bymkt(4:6,1));...
                      sum(data_sahw_t_lb(7:9,:,:).*women_moments_bymkt(7:9,1))];

data_sahw_t_byeduct_ub = [sum(data_sahw_t_ub(1:3,:,:).*women_moments_bymkt(1:3,1)); ...
                       sum(data_sahw_t_ub(4:6,:,:).*women_moments_bymkt(4:6,1));...
                      sum(data_sahw_t_ub(7:9,:,:).*women_moments_bymkt(7:9,1))];

%--- Build figure

X=[1 2 3 4 ];
hh_age_intervals = {'\leq 28', '[29,34]', '[35,40]','\geq 41'};

Y=[1 2 3 4 5 6 7 8 9];
couple_types = {'(hs,hs)', '(hs,sc)', '(hs,c+)','(sc,hs)', '(sc,sc)', '(sc,c+)','(c+,hs)', '(c+,sc)', '(c+,c+)'};

t=tiledlayout(2,2);

% hs
nexttile
plot(X,model_sahw_t_byeduct(1,:),'--', 'MarkerFaceColor',[0.5,0.5,0.5],'Color', [0 0.3470 0.5410],'LineWidth',3)
hold on
plot(X,data_sahw_t_byeduct(1,:),'-', 'MarkerFaceColor',[0.5,0.5,0.5],'Color', [0 0.3470 0.5410],'LineWidth',3)
hold on
plot(X,data_sahw_t_byeduct_lb(1,:),'--', 'MarkerFaceColor',[0.5,0.5,0.5],'Color', [0.5,0.5,0.5],'LineWidth',2)
hold on
plot(X,data_sahw_t_byeduct_ub(1,:),'--', 'MarkerFaceColor',[0.5,0.5,0.5],'Color', [0.5,0.5,0.5],'LineWidth',2)
hold on
xticks(1:9)
set(gca,'xticklabel',hh_age_intervals','FontSize',16)
xaxisproperties.TickLabelInterpreter = 'latex'; % latex for x-axis
xlabel('Fraction hs wives staying at home by household age interval','FontSize',18,'FontWeight','bold') 
ylim([0.06 0.35])

% sc
nexttile
plot(X,model_sahw_t_byeduct(2,:),'--', 'MarkerFaceColor',[0.5,0.5,0.5],'Color', [0 0.3470 0.5410],'LineWidth',3)
hold on
plot(X,data_sahw_t_byeduct(2,:),'-', 'MarkerFaceColor',[0.5,0.5,0.5],'Color', [0 0.3470 0.5410],'LineWidth',3)
hold on
plot(X,data_sahw_t_byeduct_lb(2,:),'--', 'MarkerFaceColor',[0.5,0.5,0.5],'Color', [0.5,0.5,0.5],'LineWidth',2)
hold on
plot(X,data_sahw_t_byeduct_ub(2,:),'--', 'MarkerFaceColor',[0.5,0.5,0.5],'Color', [0.5,0.5,0.5],'LineWidth',2)
hold on
legend('Model','Data', '95% CI', 'Location','northeAst','NumColumns', 3,'FontSize',14, 'FontWeight','bold') 
xticks(1:9)
set(gca,'xticklabel',hh_age_intervals.','FontSize',16)
xlabel('Fraction sc wives staying at home by household age interval','FontSize',18,'FontWeight','bold') 
ylim([0.06 0.35])

%cp
nexttile
plot(X,model_sahw_t_byeduct(3,:),'--', 'MarkerFaceColor',[0.5,0.5,0.5],'Color', [0 0.3470 0.5410],'LineWidth',4)
hold on
plot(X,data_sahw_t_byeduct(3,:),'-', 'MarkerFaceColor',[0.5,0.5,0.5],'Color', [0 0.3470 0.5410],'LineWidth',4)
hold on
plot(X,data_sahw_t_byeduct_lb(3,:),'--', 'MarkerFaceColor',[0.5,0.5,0.5],'Color', [0.5,0.5,0.5],'LineWidth',2)
hold on
plot(X,data_sahw_t_byeduct_ub(3,:),'--', 'MarkerFaceColor',[0.5,0.5,0.5],'Color', [0.5,0.5,0.5],'LineWidth',2)
hold on
xticks(1:9)
set(gca,'xticklabel',hh_age_intervals.','FontSize',16)
xlabel('Fraction c+ wives staying at home by household age interval','FontSize',18,'FontWeight','bold') 
ylim([0.06 0.35])

% duration
nexttile
plot(Y,duration_all,'--', 'MarkerFaceColor',[0.5,0.5,0.5],'Color', [0 0.3470 0.5410],'LineWidth',4)
hold on
plot(Y,duration_all_data,'-', 'MarkerFaceColor',[0.5,0.5,0.5],'Color', [0 0.3470 0.5410],'LineWidth',4)
hold on
plot(Y,duration_lb,'--', 'MarkerFaceColor',[0.5,0.5,0.5],'Color', [0.5,0.5,0.5],'LineWidth',2)
hold on
plot(Y,duration_ub,'--', 'MarkerFaceColor',[0.5,0.5,0.5],'Color', [0.5,0.5,0.5],'LineWidth',2)
hold on
xticks(1:9)
set(gca,'xticklabel',couple_types.','FontSize',16)
xlabel('Duration of marriage, by couple type','FontSize',18,'FontWeight','bold') 
ylim([1 10])

%--- Store
set(gcf,'position',[x0,y0,width,height]);
saveas(gcf, 'Output/Figure_A1.png')  
close

